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This paper begins with a general intro- 
duction to the symmetric level-index, 
SLI, system of number representation 
and arithmetic. This system provides a 
robust frameworlc in which experimen- 
tal computation can be performed with- 
out the rislc of failure due to overflow/ 
underflow or to poor scaling of the 
original problem. There follows a brief 
summary of some existing computa- 
tional experience with this system to 
illustrate its strengths in numerical, 
graphical and parallel computational 
settings. An example of the use of SLI 
arithmetic to overcome graphics failure 
in the modeling of a turbulent combus- 
tion problem is presented. The main 



thrust of this paper is to introduce the 
idea of SLI-linear least squares data fit- 
ting. The use of generalized logarithm 
and exponential functions is seen to 
offer significant improvement over the 
more conventional linear regression 
tools for fitting data from a compound 
exponential decay such as the decay of 
radioactive materials. 
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1. Introduction 

In the field of scientific computation generally 
and especially in experimental computing which is 
often at the heart of simulation and modeling prob- 
lems, the availability of a robust system of arith- 
metic offers many advantages. Many such 
computational problems are prone to failure due to 
overflow or underflow or to a lack of advance 
knowledge of a suitable scaling for the problem. 
The use of a computer arithmetic system which is 
free of these drawbacks would clearly alleviate any 
such difficulties. 

One such arithmetic is the symmetric level index, 
SLI, system. (See [3] for an introductory summary.) 
This system was developed from the original level- 
index system of Clenshaw and Olver [1] and has 
been studied in a number of subsequent papers. 
Working to any finite precision, it is closed under 



the four basic arithmetic operations (apart from di- 
vision by zero, of course) and is therefore free of 
underflow and overflow. The arithmetic system al- 
lows very large or very small numbers which may 
not be representable in a conventional floating- 
point system to be used during interim computa- 
tion while still returning meaningful results. 

Section 2 of the paper describes, briefly, the SLI 
representation and its basic algorithms and proper- 
ties. The natural error measure for SLI arithmetic, 
generalized precision, is also introduced. This mea- 
sure will be used in the SLI-linear least squares 
data fitting experiments described in Sec. 4. 

In Sec. 3, a brief summary of some of the existing 
computational evidence supporting the use of SLI 
arithmetic is presented. This concentrates first on 
two examples (the robust computation of binomial 
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probabilities and the solution of polynomial equa- 
tions by the root-squaring technique) which 
demonstrate the ability of the system to recover 
valuable information from interim results which 
would exceed the range of any floating-point sys- 
tem. The other principal example reviewed here is 
taken from the modeling of turbulent combustion. 
In this case, it is seen that commercial contour 
plotting packages used with the floating-point sys- 
tem produce seriously misleading (but initially 
plausible) results while the true situation is re- 
vealed by the use of SLI arithmetic. In this case 
there is no straightforward rescaling of the original 
problem which would allow successful floating- 
point computation. 

In Sec. 4, we concentrate on a topic which is 
commonly used in simulation and statistical 
analysis— least squares data-fitting. Our interest is 
in curve-fitting, not the related problem of parame- 
ter estimation. The particular problem discussed is 
a compound exponential decay such as might be 
encountered in modeling the decay of radioactive 
materials. Compound decays with widely varying 
half-lives are not well approximated by the com- 
monly used log-linear least squares method. The 
use of the SLI representation function— a general- 
ized logarithm— permits much better agreement 
with the data while still using very low degree ap- 
proximating functions. The benefits derived from 
using SLI arithmetic— or, in this case, just the SLI 
representation— are made plain by a sequence of 
graphical examples. 

2. The SLI Arithmetic System 

The level-index number system for computer 
arithmetic was first suggested in Clenshaw and 
Olver [1], [2]. The scheme was extended to the 
symmetric level index, SLI, representation in [4] 
and has been studied in several further papers in 
the last few years. Much of the earlier work is sum- 
marized in the introductory survey [3]. The primary 
virtue of SLI arithmetic is its freedom from over- 
flow and underflow and the consequent ease of al- 
gorithm development available to the scientific 
software designer. This is not the only arithmetic 
system that has been proposed with this aim: for 
example, the work of Matsui and Iri [14], Hamada 
[7], [8], and Yokoo [27] suggested and studied 
modified floating-point systems which share some 
properties of level-index; Smith, Olver, and Lozier 
[21] studied an extended range arithmetic. 

Possible hardware implementations of SLI arith- 
metic were discussed in [18], [23], and [26] while a 



software implementation incorporating some ex- 
tended arithmetic was described in [25]. The error 
analysis of SLI arithmetic is discussed in [2] and [4] 
and is extended in [11], [16], and [17]. Applications 
and software engineering aspects of the level-index 
system have been discussed in [5], [10], [12], and 
[24]. 

The SLI representation of a real number X is 
given by 

where the two signs ^^and rx are ± 1 and the gener- 
alized exponential function is defined for j; > by 



^^ ^ \exp{^x-l)) 



x>l 



It follows that for Z>1, 



^ = exp(exp(. ..(exp/)...)) 

where the exponentiation is performed l=[x] times 
and x=l+f. The integer part / of a: is called the 
level and the fractional part / is called the index. 

The freedom of this system from over- and un- 
derflow results from the fact that, working to a pre- 
cision of no more than 5 million binary places in 
the index, the system is closed under the four basic 
arithmetic operations with only three bits allotted 
to the level. This is discussed briefly in [1], [4] and 
considered in some detail in [11]. 

The basic SLI arithmetic operation is that of 
finding the SLI representation Sz<f>(zy' of Z =X ±Y 
where X, Y are also given by their SLI representa- 
tions. Without loss of generality, we may assume 
that X^Y>0 so that 5j = +l. The computation 
entails the calculation of the members of three 
short sequences which vary according to the partic- 
ular circumstances. In every case, the sequence de- 
fined by 



a,= 



1 



' Hx-j) 



(J =1-1, 1 -2 0) 



where l = [x], is computed using the recurrence 
relation 



_--f 



ay-i = exp(-l/ay); a/-i = e 

Depending on the values of r^, ry, and Tj, the other 
sequences that may be required are given, for ap- 
propriate starting values, by 
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"'-^ix-jy'^ <l>(y-jy "> <l>(y-j) 



^'=^' hj = Hz-J). 



Specifically, the sequence (bj) is used when 
fy = + 1. The terms can be computed using 

fe;-. = exp(^) 

with the initial value given by 

u - ^_ fexpfe-1/flm) (m<l) 
ftm-i-flm-ie^- lexpfe-/) (m=/) 

where m =[y] is the level ofy. Since, in this case, 
y^x, it follows that 0^bj<l. 

The sequence (uj) is used when rx = +l,ry=—l 
and is computed like (fly). It is similarly bounded: 
O^fl;, aj<l. 

For the case where rx=ry = —l, the requirement 
^>y implies ;c^}'. The sequence (oj) is computed 
as before along with the sequence (/3;). This latter 
is computed using the recurrence relation 

for J </ with the initial value 

a _fexp(f-l/a/) (Km) 
P'-'-Uxp(f-g) ^i-m)- 

The results of these calculations can be combined 
to yield, for the first two cases 

Co = l±feo or Co=l±aoao 



while for the "small" case, we compute 



co'=-=l±i3b 

Co 



from which the required terms of the c -sequence 
may be computed by the recurrence 

Cj =1+0; In C;-l . 



The number of such terms is bounded by / — 1 after 
which some terms of the sequence hj may be 
needed. This just involves forming repeated natural 
logarithms of c/-i. Detailed derivations of the vari- 
ous sequences can be found in [2] and [4]. 

This algorithm is implemented in the software 
implementation of SLI arithmetic described in [25]. 
This implementation also takes advantage of the 
relative ease of performing extended arithmetic op- 
erations such as summation or forming scalar prod- 
ucts in SLI arithmetic, the details of which are 
discussed more fully in [25], [26]. The major advan- 
tage of these extended operations is likely to be in 
parallel computing environments where most of 
the likely speed loss suffered by SLI arithmetic will 
be recouped. 

The results of [26] indicate that a parallel SLI 
processor could yield a reasonable speed-up over 
serial floating-point computation for extended 
operations of even moderate length. A fairer "par- 
allel-parallel" comparison suggests a likely slow- 
down of arithmetic by a factor of around 2 for 
extended operations— which probably represents a 
loss of some 10 to 20% in run-time. However, even 
that small price to pay for a robust arithmetic is 
likely to be recouped for many parallel operations. 

In [12], it is seen that several basic operations 
such as forming vector norms are difficult to pro- 
gram both efficiently and robustly for floating- 
point machines but become straightforward tasks 
in the SLI environment. It is the simplicity of pro- 
gramming, and code which is free of special scaling 
or exception-handling cases, which are likely to tip 
the balance in favor of SLI arithmetic once suitable 
hardware implementations exist. 

The appropriate error measure for computation 
in the level index system is no longer relative error 
(which corresponds approximately to absolute pre- 
cision in the mantissa of floating-point numbers) 
hvA generalized precision which corresponds to abso- 
lute precision in the index. This error measure is 
introduced in [1]. Generalized precision has some 
significant advantages compared to relative error, 
not least of which is that it is a metric so that the 
symmetry of jc approximating jc, and jc approximat- 
ing ;c, is a natural aspect of the error analysis. De- 
tailed error analyses of numerical processes will 
inevitably be different in this system than for any of 
the floating-point systems but significant progress 
has already been made in this respect. (See, for 
example, [1], [17], [26].) 
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Again considerable benefits can be achieved for 
extended calculations. Olver [17] has demonstrated 
that, at relatively low cost, it is possible to perform 
a concurrent error analysis. Such analysis is partic- 
ularly well-suited to a parallel environment since it 
would be performed by simultaneous duplication of 
the operations for slightly adjusted data. The ad- 
justments are similar to the use of directed round- 
ing in interval arithmetic and have a similar effect 
in yielding guaranteed error bounds for the results 
obtained. 

A first-order error analysis for extended sums 
and products [26] leads to conclusions that are 
broadly similar to those for the floating-point sys- 
tem. However, for the SLI system, we find that the 
generalized precision of the final result is bounded 
by Nil times the generalized precision for individ- 
ual operations. 

3. Computational Evidence 

We give three examples of computations which 
are highly susceptible to underflow and/or over- 
flow, and we present computational evidence that 
SLI arithmetic produces valid results. The first ex- 
ample, the binomial probability distribution, was 
treated in [22] to introduce and evaluate tech- 
niques for dealing with underflow and overflow in 
floating-point, FLP, computation. The second ex- 
ample, the classical Graeffe algorithm for deter- 
mining the zeros of a polynomial, is almost never 
used in FLP computation because the root-squar- 
ing process almost always introduces very large 
and/or very small numbers as intermediate results. 
The final example is drawn from a computer 
graphics display of an analytical solution to a 
model problem in turbulent combustion theory. 

3.1 Binomial Probability Distribution (BPD) 

In addition to [22] this example was treated in 
[14] and [3]. The BPD and its importance in statis- 
tics is discussed, for example, in [15], pp. 54-58. 
Let /J be the probability of a favored outcome from 
an individual chance event. Then q = \-p is the 
probability that the favored outcome will not occur. 
Now, the probability of the favored outcome occur- 
ring exactly k times in n events is 



/-a) 



/^"-* 



kc=[{n+\)p\ 

and/t decreases steadily as k moves away from kc. 
The values of/* are exceedingly small when n is 
of moderate size and k is not close to kc. For exam- 
ple, [22] considers n =2000 and/j =0.1. Since the 
underflow threshold' in IEEE standard 32-bit FLP 
arithmetic is 1.18x10"^^ for normalized numbers, 
fk underflows when 0«/:«51 and 393«/:<2000. 
Replacement of these underflows by zero (or de- 
normalized numbers) may be acceptable for many 
purposes, such as for the computation of the cumu- 
lative probability 

m 

I{n,m,p)=^fk (O^m^n). 

In reality, however,/* and / are always strictly posi- 
tive, and it can be disconcerting when a computer 
program prints zero. 

A more serious problem, perhaps, is that an FLP 
algorithm must guard carefully against intermedi- 
ate underflow and overflow. Overflow is possible in 
forming the binomial coefficient 



(nyjn 



-l),..(n-k + l) 



k\ 



and underflow is possible in forming powers of p 
and q. Indeed, with n =2000, p =0.1 and k =k, = 
200, 

/2oo= 0.02972 287717 
but 

[l) = 10^"-^, p^'^= 10-^, q"^ = lO-''^-^'^ 

Since </* < 1, the algorithm cannot be allowed to 
fail because of overflow. Zero is an acceptable re- 
sult only if fk is below the underflow threshold. 

Two algorithms are considered in [22]. Al- 
gorithm I forms /i by (i) multiplying in all factors of 
the numerator of the binomial coefficient, (ii) di- 
viding out all factors of the denominator, (iii) mul- 
tiplying in all factors of /;*, and (iv) multiplying in 
all factors of ^""*. As we have seen, this algorithm 
is sure to fail in many cases of interest. It is ren- 
dered usable in [22] by introducing a counter / and 
two large positive constants ci and cz such that ciCa 
is slightly less than the overflow threshold and 1/ 
CiCz is slightly more than the underflow threshold. 



This distribution attains its maximum when k 
equals its modal value 



' Here and elsewhere in this section, numbers written with a 
decimal point are correct to all places shown. 
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Let us refer to this modification as Algorithm lA. 
Before each arithmetic operation, Algorithm lA 
tests the current result. If it does not lie between 
1/ci and ci, it is scaled into this interval by multiply- 
ing or dividing by ci and incrementing or decre- 
menting / accordingly. At the end, cifk is represent- 
ed as a normalized FLP number with / 5 0. If / = 0, 
the algorithm is complete. If / = ! and ci/ji<l/c2, 
the algorithm divides once by ci and returns fk . In 
all other cases the algorithm returns zero. 

Algorithm II forms fk from the same operations 
as Algorithm I but in different order. Operations 
from part (i) increase the current result, whereas 
those from the other three parts decrease it. Al- 
gorithm II uses only a single large constant cu It 
starts with increasing operations. When the current 
result would exceed Ci, Algorithm II switches to de- 
creasing operations until the current result would 
fall below 1/ci, at which point it switches back to 
increasing operations. The algorithm ends when ei- 
ther fk is completely formed or the current result 
falls below 1/ci and no increasing operations re- 
main, in which case fk is returned as zero. 

These algorithms may be compared as follows. 
Algorithm I produces the shortest program with 
the simplest logic. It requires no tests or system-de- 
pendent constants. It is feasible for SLI but not for 
FLP arithmetic. Algorithm lA produces a program 
more than twice as long because, at every arith- 
metic operation, it requires a test and contingent 
coding to handle underflow or overflow. The num- 
ber of contingent operations depends onp,n,k,Ci 
and cz. The latter two constants depend on the un- 
derflow and overflow thresholds but there is no 
natural mathematical definition of this depen- 
dence. A counter is required also. The algorithm is 
designed for FLP and is not appropriate for SLI 
arithmetic. Algorithm II produces a program that 
is slightly more complex and slightly shorter than 
Algorithm lA. It requires a test but no contingent 
coding at every arithmetic operation, only one con- 
stant, and no counter. It is feasible for both FLP 
and SLI arithmetic. 

Table 1 summarizes these observations. The 
table also gives some data on the relative error in 
the computed value of /200, with p =0.1 and 
n =2000, in 32-bit SLI and FLP arithmetic. Mea- 
sured relative errors were obtained by comparison 
against 64-bit FLP calculations. The measured er- 
ror is similar for Algorithm II in both arithmetics, 
with the smaller error occurring in SLI. For Al- 
gorithm I, the FLP error is » because of overflow 
failure. The SLI error is larger than for Algorithm 
II because the simpler algorithm generates larger 



intermediate values (up to 10*^'-^ for Algorithm I, 
10^*^ for Algorithm II). Algorithm lA is not appro- 
priate for SLI and was used only as a specific rem- 
edy for FLP arithmetic. It performed 31 contingent 
operations in computing /200 and produced an error 
slightly larger than the FLP error of Algorithm II. 



Table 1. Attributes of three algorithms for the binomial proba- 
bility distribution. It = number of events, /:= number of favored 
outcomes, C = number of contingent operations to avoid under- 
flow and overflow 





I 


lA 


II 


Program length 


Short 


Long 


Long 


Program logic 


simple 


Complex 


Complex 


No. constants 





2 


1 


No. counters 





1 





No. operations 


n+2k 


n+2k+c 


n+2k 


No. tests 





n+2k 


n+2k 


Measured rel err 


6.4X10-" 


N.A. 


-1.5x10-5 


(SLI) 








Measured rel err 


00 


-4.5x10-' 


-4.3X10-' 


(FLP) 









3.2 The Graeffe Root-Squaring Process 

At the end of the preceding subsection, it was 
suggested that one algorithm for the BPD leads to 
larger relative errors than another because the first 
algorithm generates larger intermediate values. Al- 
though this is true for the BPD, it is not a reliable 
guide for all algorithms. Indeed, it is easily proved 
that the relative error in a function _y =/(ji:) caused 
by a relative error its argument is approximated to 
first order by 



where 



8y='Af{x)8x 



Afix) = 



- xf'jx) 



fix) 



This function is sometimes called the relative-error 
amplification factor for /. However, it either ampli- 
fies or deamplifies the relative error according to 
whether \Af(x)\ > 1 or \Af(x)\ < 1. If \(Afix)\<h the 
deamplification effect is very strong. An example is 
y ^x" for 0<a <L Here the deamplification factor 
isyl/(ji;)=a. 

The Graefi'e root-squaring process is a very old 
numerical method for solving algebraic equations; 
in [9] it is traced back to 1762. For more modern 
accounts, see [6], [9], and [20], pp. 1174-8. Let 

p{x)=aX ^an-^x"'"^ ^ . . . +ao, a«5^0 
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be an arbitrary polynomial. If y is a zero of p{x), 
the Graeffe process forms approximations to lyl"* 
where m can be arbitrarily large. More specifically, 
if a subset of zeros lies on or very near a circle in 
the complex plane, the method finds the mth 
power of the radius of the circle and an algebraic 
equation for the mth powers of the zeros in the 
subset. The size of m in any particular application 
depends on how well separated these circles are - 
the smaller the separation, the higher m must be. 
The idea of the method is to increase the separa- 
tion of the circles by squaring their radii. 

The Graeffe process involves the construction 
[20], p. 1174 of a finite number of the polynomials 
in the sequence 

p<-''\x)=aPx"+a^'lix"-' + ...+aV,a'^?^0 

such that the zeros of/>'''(x;) are the zeros ofp(x) 
raised to the 2''th power. Here />'°'=/>— that is, 
a '■f = aj for all ; — and 

min (j,n -j) 

«/+"=(- iy[flWp+2 2 (-i)V-w/^*. 

In particular. 

As the iteration proceeds, some (or none or all) of 
the a '/'with 0<j<n begin to satisfy the approxi- 
mation 

a.(r+i)=(_iy[fl(j-)]2. 

Let us suppose this happens for some r and 

0<;i<;2 <...<;.<« 

and let us define ;o = 0, ;j+i=n. Then there are 
zeros oi p(x) that lie on or very near a circle of 
radius 

lr.| = ki?/«>iir(i = i,2,...,s+i). 

Because the deamplification factor 2"' is small, this 
determination of |y,| should be very stable. This 
expectation is borne out by examples in [6], [9], and 
[20], all done by desk calculator or slide rule in the 
era before FLP computation became widespread. 
Indeed, on p. 187 of [6] it is stated that Graeffe's 
method is "generally useful" and "well adapted to 
the computing machine or slide rule"; obviously 
the author did not appreciate the impact of under- 
flow or overflow! In most of these examples num- 
bers are generated that significantly exceed the 



IEEE overflow threshold, even in double precision. 

When,y =n -1, the preceding discussion suffices 
to summarize the use of the Graeffe process up to 
the determination of the phase angle; if the polyno- 
mial has all real coefficients, the phase angle is ei- 
ther or TT and the determination can be made by 
substitution of ± y, into the original equation. In all 
other cases the method needs modification or ex- 
tension. Some of these are indicated in [6], [9], and 
[20], including the important case of real polynomi- 
als with nonrepeated real and complex conjugate 
zeros. 

The Graeffe process has fallen into disuse be- 
cause it is not well suited to FLP computation. Cer- 
tainly alternative methods are available to solve 
algebraic equations. One popular method is to use 
an algorithm from numerical linear algebra to find 
the eigenvalues of the so-called companion matrix 
of p(x). However, the Graeffe process could 
equally well be applied to the problem of finding 
matrix eigenvalues. This is a powerful motivation 
for considering SLI computer arithmetic. 

Such a consideration was begun in [3] and [5], 
where the method is described and a preliminary 
error analysis is given. Several more substantial ex- 
amples are computed in SLI arithmetic with very 
satisfactory results even when very large values of 
m were needed. The method has also received at- 
tention in [14]. 

3.3 Graphics for a Combustion Problem 

A model problem in the theory of turbulent com- 
bustion, introduced in [13], involves two chemical 
species (fuel and oxidizer) which occupy two adja- 
cent half-spaces separated conceptually by an im- 
pervious plane boundary. At time t=Q, the 
boundary loses its imperviousness and a line vortex 
is imposed in the plane of the boundary, leaving 
the two species free to diffuse into each other, re- 
act, and produce a flame surface at the boundary, 
which is distorted by the vortex into a cylindrically 
symmetric spiral surface. With the axis of the vor- 
tex identified as the z-axis, cylindrical symmetry re- 
duces the spatial variables to plane polar 
coordinates. The location and shape of the flame 
surface is of particular interest in applications. For 
example, reactant consumption and heat produc- 
tion are obtained by integrating the normal deriva- 
tive of the solution along the flame surface. 

This problem has undergone extensive analytical 
and computational development in recent years. 
For example, the development in [19] of a solution 
in terms of a similarity variable 
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^ = te 



where v is the kinematic viscosity, opened up the 
possibility of a two-dimensional analysis. 

Indeed, an explicit representation was obtained 
of a function Z (tj , ©) as a Fourier series. The flame 
surface of the model problem corresponds to the 
level curves, or contours, of this function. There- 
fore, it can be computed (in principle, at least) by 
inverse interpolation from data on a grid. The data 
are the computed values of Z on the grid. Powerful 
graphics software is available to compute and dis- 
play contours on a wide variety of graphics devices. 
Unfortunately, this software is not robust in the 
face of underflow, as we shall see in the case of the 
flame surface. 

From a mathematical as well as a scientific 
standpoint, Z(ri,6) exhibits interesting behavior. 
As 77 -*0 for fixed 6, it oscillates with unboundedly 
growing frequency while its magnitude tends to 
zero exponentially fast. Of particular interest is the 
shape of the contours Z= ±e. For a small e>0, 
the two contours reflect each other in the origin. 
Their shape is a spiral which winds toward the 
origin up to a point that depends on e, at which it 
abruptly reverses direction and winds away from 
the origin. In the limit as e-*0, the contours fuse 
into one that can be regarded as having two 
branches, one a spiral that winds right down to the 
origin while encircling it infinitely often and the 
other its reflection in the origin. 

This complicated behavior suggests that it would 
be a difficult challenge for graphics software to 
produce correct contours. Indeed, this is the case. 
IEEE single-precision arithmetic, with a word 
length of 32 bits, was used to execute a standard 
library subroutine for contour plotting. Fig. 1 shows 
the results for e = 0. Clearly, the central contour 
does not exhibit the required symmetry or spiral 
behavior. 

Figure 2 was produced in exactly the same way 
except SLI arithmetic was used. In particular, no 
algorithms were changed and no scaling was intro- 
duced. The mesh of 151^ points is identical. The 
smallest magnitude of data on the mesh, excluding 
the origin, was 10"''^'^, approximately. The contour 
shows the correct qualitative and quantitative be- 
havior within the limits of the mesh resolution. 
With a refinement of the mesh, the infinite spiral 
would be resolved correspondingly closer to the 
origin; in fact, this was demonstrated on a mesh 
nine times finer, in which case the smallest magni- 
tude was lO"^****^, approximately. 




Fig. 1. Contours of 2 = - 0.4(0.1)0.4 computed in IEEE single 
precision by the GCONTR grapliics subroutine. Resolution = 

15P. 




Fig. 2. Contours of 2 = - 0.4(0.1)0.4 computed in 32-bit SLI by 
the GCONTR graphics subroutine. Resolution = 151^. 

Figure 3 is a repeat of Fig. 1 but with one modi- 
fication: A black mark is plotted at every mesh 
point where the functional data was zero. The only 
mesh point where Z = is the origin; all the other 
zeros are the result of underflow. This shows clearly 
that the sole cause of the graphics failure in FLP 
arithmetic is underflow. 
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Fig. 3. Contours of Z = - 0.4(0.1)0.4 computed in IEEE single 
precision by the GCONTR graphics subroutine. Resolution = 
151^. Blacic marks indicate points where Z underflowcd. 



From a scientist's point of view, the meaning of 
these results is the following. Small values of c cor- 
respond to cases of near stoichiometry among the 
concentrations of the combustion reactants. Perfect 
stoichiometry corresponds to complete combustion, 
in which all reactants are consumed fully. If a nu- 
merical experiment is conducted in which e takes 
the sequence of values 10"^ 10"*, 0, say, the results 
are confusing. That is, the scientist sees correct be- 
havior for the nonzero values but incorrect behavior 
for e = 0. The transition to perfect stoichiometry 
cannot be made visible graphically with any cur- 
rently available floating-point processor. This de- 
fect does not arise with SLI arithmetic. 

We conclude this section by remarking that scal- 
ing does not provide an effective solution to the 
FLP graphics failure. The stoichiometric contour 
Z = is obtained by inverse interpolation from 
nearby data on the grid. For grid points near the 
origin, as we have seen, this data is exceedingly 
small. However, for grid points near the boundary 
of the plotting region, the data is not small; cer- 
tainly it is not below the underflow threshold. The 
data varies over some 4O0O0 orders of magnitude. 
Since IEEE single precision varies over only 76 or- 
ders of magnitude, approximately 526 scalings 
would be required. Clearly, scaling attacks the fail- 
ure only very weakly. Added to this weakness would 
be the difficulties presented by the contouring soft- 



ware itself. The scaling in this example needs to be 
done in thin annuli around the origin but the soft- 
ware works with rectangular regions. Therefore, a 
complicated mosaic of subregions would need to be 
combined to form the complete contour plot. 

4. SLI-Linear Least Squares Data-Fitting 

This section is concerned with a common tech- 
nique in simulation — least squares fitting of data. 
The experiments reported provide a comparison 
between linear least squares, log-linear least 
squares and SLI-linear least squares as curve fitting 
techniques for data from compound exponential de- 
cay functions such as might arise in studying ra- 
dioactive decay. 

The principle of linear least squares polynomial 
fitting to data is to find that polynomial of fixed 
(maximum) degree, N say, for which the sum of 
squared errors 

n 
1=0 

is minimized. Here, the data consists of the ordered 
pairs (Xi,yi) and it is assumed that the number of 
data points n+1 exceeds A'^ (usually by a significant 
proportion). 

The log-linear fit to this data uses the approxima- 
tion exp(qN(x)) where qN is the least squares poly- 
nomial approximation to the data (Xi,\nyi). 

The SLI-linear fit uses <P{qff{x)) where, this 
time, qu is the least squares polynomial approxima- 
tion to the data (x, , '^(yi )) where the functions (?, ^ 
are the functions used in the SLI representation of 
positive real numbers. They are defined by: 



nx)= 



r iifix)-i x^i 



and 



^(.Uil+x) x^O 
^^""^ U(l + |;c|)-' x<0 



where ^ is the generalized exponential defined in 
Sec. 2 and the generalized logarithm \f/ is its inverse. 
The first experiments all used the same test func- 
tion to generate the data on [0, 5]: 

/(x) = 20e-^^ + 5e-" 
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which consists of two exponential decay compo- 
nents the first of which is initially much the greater 
but decays much faster. In all cases, the least 
squares approximations of the appropriate degree 
using polynomial, log-linear and SLI-linear fitting 
are plotted together with the original function for 
comparison. 

The first set of graphs^ Fig. 4 uses (degree 1) 
linear least squares fits to randomly generated data 
sets in [0, 5] subject only to a bias towards zero in 
the selection. This was achieved by forcing approxi- 



mately 50% of the data to lie in [0, 1], 25% in [1, 
2], 12.5% in [2, 3] and so on. For the first pair of 
graphs the choice of data points was entirely ran- 
dom save for the distribution described above. The 
best fit was achieved with the SLI-linear approxi- 
mation. In the second pair. Fig. 5, the further re- 
striction that jcq^O was enforced. Again the 
SLI-linear approximation is the best in both cases 
and forcing to be a data point has improved sig- 
nificantly the "goodness of fit" near the origin. 





(a) 



(b) 



Fig. 4. (a) 8 data points, and (b) 32 data points. 





(a) 



(b) 



Fig. 5. (a) 8 data points including 0, and (b) 32 data points including 0. 
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The next set of graphs, Fig. 6, shows the result of 
quadratic approximations using diiferent numbers 
of randomly chosen data points with the constraint 
that two of them were tied to the end-points. In 
every Instance two of the curves were almost indis- 
tinguishable: the test function and the SLI fit. The 
polynomial fit was uniformly the worst— as would 
be e^qjected. Typically, the logarithmic fit either re- 
mains too high throughout the region or starts too 
low and then crosses to become too large. The sec- 



ond parr of graphs in Fig. 6 magnifies the region 
[0,1.5] X [0,15] for different random data sets. 

The effect of an entirely random data set is illus- 
trated in Fig. 7. With even three data points the 
SLI-fit again produces tolerably good agreement 
with the original function— in sharp contrast to 
both the others having reached their respective 
minima to the left of x =2. With a larger number of 
data points, the SLI approxnnation again repro- 
duces the original function to high accuracy. 





(a) 



(b) 





(c) 



(d) 



Fig. 6. (a) 3 points, (b) 40 points, (c) 3 points, magnified view, and (d) 40 points, magnified view. 
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(a) 



(b) 



Fig. 7. (a) 3 points, and (b) 35 points. 



The observation in the Figs. 6 and 7 that the 
SLI-quadratic fit using just three points (including 
the end-points) has produced such close agreement 
with the test function suggested further investiga- 
tion to see whether this low degree approximation 
can be relied upon to produce high-accuracy ap- 
proximations. This reliability persists for all choices 
of three data points {0, ill, 5} for / = 1, 2, . . ., 9. 
The results for i — l and 9 are presented in Fig. 8. 
Again the second pair magnifies the region 
[0,1.5] X [0,15]. 

Consideration of a table of differences suggests 
that interpolation at just two points may yield sur- 
prisingly good results. The evidence thus far is that 
use of the left-hand end-point is advantageous. 
This proved to be especially true in the linear inter- 
polation case. Interpolation at and / for / = 1, . . . , 
5 continued the pattern of the SLI fit being visibly 
superior. The graphs for the cases i—l and 4 are 
reproduced in Fig. 9. 

Experiments using 3 components including a 
term with much slower decay again had the consis- 
tent result that the SLI-approximations were supe- 
rior. See Fig. 10. The data function used was: 



/(:«:) =20e-^+5e-^ + e" 



■xIA 



The first of the graphs in Fig. 10 shows the result of 
a quadratic approximation using 20 random data 
points including the end-points. The SLI approxi- 
mation is almost exact except at the left-hand end 
where, like the others, it lies below the data curve. 
The second one demonstrates that increasing the 
degree of the approximation to cubic yields almost 
perfect agreement with a mere 7 data points. The 
polynomial fit has its local maximum around a: = 4. 
The log-linear fit is also effective beyond about 
a: = 2 by which stage the two faster decaying terms 
have become insignificant. 

The final example was to approximate a "faster- 
than-exponential" decay given by 

f{x) = 15tx^{-x'^-\ll) + 5 exp(-A:) 

over the interval [0, 2.5]. The graph is typical of the 
results which again indicate a superiority for the 
SLI-approximation. In this case quadratic approxi- 
mations and 20 randomly placed data points were 
used. The SLI approximation is the only one with 
the right "shape" having an inflection point very 
close to that of the test function. The results are il- 
lustrated in Fig. 11. 
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(a) 



(b) 





(c) 



(d) 



Fig. 8. Three fixed data points: endpoints and stated point, (a) 0.5, (b) 4.5, (c) 0.5, magnified view, and (d) 4.5, magnified view. 
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(a) 



(b) 



Fig. 9. (a) Data points 0, 2, and (b) data points 0, 4. 



25.DD 



D.DD 




D.DD 



' 1 11,-^ , — 



5.DD 




(a) 



(b) 



Fig. 10. (a) 20 points, quadratic fit, and (b) 7 points, cubic approximation. 
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Fig. II. Faster than exponential decay, 20 random data points, 
quadratic approximations. 



5. Conclusions 

The principal conclusions to be drawn from this 
paper are that SLI arithmetic offers a robust alter- 
native to the floating-point system which enables 
many of the standard tasks of scientific computing 
to be performed in a simple yet reliable way. The 
computational evidence of Sec. 3 summarizes its 
power in three inherently different situations: com- 
puting binomial probabilities, solving polynomial 
equations and contour plotting for functions which 
have widely varying values and so are not amenable 
to scaling. 

The final section discusses the use of SLI arith- 
metic in a curve-fitting situation and demonstrates 
that the SLI representation can be used to advan- 
tage in fitting data from a compound exponential 
decay. The experimental evidence indicates that a 
very good fit can be obtained with a low order SLI- 
linear least squares fit using only a small number of 
data points. This topic will be the subject of further 
experiment and analysis. 

6. References 

[1] C W. Clenshaw and F. W. J.OIver, Beyond floating point, 
J. ACM 31, 319-328 (1984). 

[2] C W. Clenshaw and F. W. J.OIver, Level-Index arithmetic 
operations, SIAM J. Num. Anal. 24, 47(M85 (1987). 

[3] C. W. Clenshaw, F. W. J. Olver, and P. R.Turner, Level- 
index arithmetic: An introductory survey. Numerical 
Analysis and Parallel Processing, P. R. Turner, ed.. Lec- 



ture Notes in Mathematics 1397, Springer Verlag (1989) 
pp. 95-168. 
[4] C. W. Clenshaw and P. R. Turner, The symmetric level- 
index system, IMA J. Num. Anal. 8, 517-526 (1988). 
[5] C. W. Clenshaw and P. R. Turner, Root-squaring using 

level-index arithmetic. Computing 43, 171-185 (1989). 
[6] N. B. Conkwright, Introduction to the Theory of Equa- 
tions, Ginn and Company (1957). 
[7] H. Hamada URR: Universal representation of real num- 
bers, New Generation Computing 1, 205-209 (1983). 
[8] H. Hamada, A new real number representation and its 
operation, Proc. ARITH8, M. J. Irwin and R. Stefanelli, 
eds., IEEE Computer Society, Washington DC, May 
1987, pp. 153-157. 
[9] C A. Hutchinson, On Graeffe's Method for the Numeri- 
cal Solution of Algebraic Equations, American Mathe- 
matical Monthly 42, 149^-161 (1935). 

[10] D. W. Lozier, An Underflow-induced Graphics Failure, 
to be published. 

[11] D. W. Lozier and F. W. J. Olver, Closure and precision in 
level-index arithmetic, SIAM J. Num. Anal. 27, 1295-1304 
(1990). 

[12] D. W. Lozier and P. R. Turner, Robust Parallel Computa- 
tion in Floating-Point and SLI Arithmetic, to be pub- 
lished in Computing. 

[13] F. E. Marble, Growth of a diffusion flame in the field of a 
vortex, in Recent Advances in Aerospace Sciences, 
C. Casci, ed., (1985). 

[14] S. Matsui and M. Iri An overflow/underflow-free float- 
ing-point representation of numbers, J. Information 
Proc. 4, 123-133 (1981). 

[15] A. F. Mood, Introduction to the Theory of Statistics, 
McGraw-Hill (1950). 

[16] F. W. J. Olver, A new approach to error arithmetic, 
SIAM J. Num. Anal. IS, 368-393 (1978). 

[17] F. W. J. Olver, Rounding errors in algebraic processes - 
in level-index arithmetic, Proc. Reliable Numerical Com- 
putation, M. G. Cox and S. Hammarling, eds., Oxford, 
(1990) pp. 197-205. 

[18] F. W. J. Olver and P. R. Turner, Implementation of level- 
index arithmetic using partial table look-up, Proc. 
ARrrH8, M. J. Irwin and R. Stefanelli, eds., IEEE Com- 
puter Society, Washington, DC (1987) pp. 144-147. 

[19] R. G. Rehm, H. R. Baum, D. W. Lozier, and J. Aronson, 
Diffusion-controlled reaction in a vortex fleld. Combus- 
tion Sci. Technol. 66, 293-317 (1989). 

[20] K. Rektorys, Survey of Applicable Mathematics, The MIT 
Press (1969). 

[21] J. M. Smith, F. W. J. Olver and D. W. Lozier, Extended- 
range arithmetic and normalized Legendre polynomials, 
ACM Trans. Math. Softw. 7, 93-105 (1981). 

[22] P. H. Sterbenz, Floating-Point Computation, Prentice- 
Hall (1974). 

[23] P. R. Turner, Towards a fast implementation of level- 
index arithmetic, Bull. IMA 22, 188-191 (1986). 

[24] P. R.Turner, Algorithms for the elementary functions in 
level-index arithmetic. Scientific Software Systems, M. G. 
Cox and J. C. Mason, eds.. Chapman and Hall (1990) pp. 
123-134. 



484 



Volume 97, Number 4, July-August 1992 

Journal of Research of the National Institute of Standards and Technology 

[25] P. R. Turner, A software implementation of SLI arith- 
metic, Proc. ARITH9, M. D. Ercegovac and E. Swartzlan- 
der, eds., IEEE Computer Society, Washington, DC, 
September 1989, pp. 18-24. 

[26] P. R. Turner, Implementation and analysis of extended 
SLI operations, Proc. ARITHIO, P. Komerup and D. W. 
Matula, eds., IEEE Computer Society, Washington, DC, 
June 1991, pp. 118-126. 

[27] H.Yokoo, Overflow/underflow-free floating-point number 
representations with self-delimiting variable length expo- 
nent field, Proc. ARITHIO, P. Kornerup and D. W. 
Matula, eds., IEEE Computer Society, Washington, DC, 
June 1991, pp. 110-117. 



About the Authors:: Daniel W. Lozier is a mathe- 
matician in the Applied and Computational Mathe- 
matics Division of the NIST Computing and Applied 
Mathematics Laboratory in Gaithersburg, MD. Peter 
R. Turner is a professor in the Mathematics Depart- 
ment of the U.S. Naval Academy in Annapolis, MD. 
The National Institute of Standards and Technology is 
an agency of the Technology Administration, U.S. 
Department of Commerce. 



485 



